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PREFACE 


The  work  we  report  here  has  two  primary  objectives.  The  first  goal  is  to  develop  methods  for 
modeling  regional  wave  propagation  in  laterally  varying  media,  and  the  second  is  to  develop 
numerical  and  theoretical  models  for  radiation  of  elastic  waves  from  explosion  sources.  Our 
report  includes  a  preprint  entitled  ’’Seismo-acoustic  multiple  scattering;  A  comparison  of 
multiple  multipole  expansions  with  an  ultrasonic  experiment,”  which  has  been  submitted 
to  the  Journal  of  the  Acoustical  Society  of  America.  This  work  presents  a  calibration  of 
the  multiple  multipole  technique  with  laboratory  data  that  we  collected,  and  demonstrates 
the  validity  of  the  approach.  Potential  applications  include  the  modeling  of  path  scattering 
effects  on  regional  seismograms  and  the  effects  of  heterogeneity  in  the  vicinity  of  the  source. 
The  second  section  of  this  report  is  a  preprint  of  the  paper  ’’Seismic  radiation  from  explosively 
loaded  cavities  in  isotropic  and  transversely  isotropic  media,”  accepted  for  publication  in  the 
Bulletin  of  the  Seismological  Society  of  America.  It  demonstrates  that  small,  decoupled 
explosions  in  cylindrical  cavities  will  generate  shear  waves  with  amplitude  depending  on 
cavity  shape  (ratio  of  length  to  width).  In  addition,  we  present  the  interesting  result  that 
locating  such  a  cavity  in  transversely  isotropic  media  can,  in  special  cases,  actually  produce 
less  shear  wave  energy  than  in  an  isotopic  medium. 
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Abstract 


The  wavefield  scattered  from  only  three  solid  rods  submerged  in  water  is  very  complex. 
Although  the  major  events  can  readily  be  interpreted  as  direct  reflections,  the  remaining 
wavefield  consists  of  complicated  interactions  of  the  wavefields  with  the  scatterers  and  the 
water  surface.  Forward  modelling  of  the  situation  is  a  viable  interpretative  aid,  especially 
if  the  method  decomposes  the  wavefields  by  their  respective  origins.  Multiple  multipole 
expansions  have  already  been  used  as  a  versatile  tool  to  model  acoustic  or  elastic  multi¬ 
ple  scattering  problems  where  homogeneous  scatterers  were  embedded  in  a  homogeneous 
fuUspace.  In  the  present  paper,  the  scheme  is  expanded  to  multiple  scattering  between  solids 
submerged  in  a  fluid.  For  each  scatterer,  the  waves  induced  in  the  fluid  are  expressed  by  sets 
of  multipole  solutions  with  different  origins.  Thus  by  construction,  the  scattered  wavefields 
are  decomposed  by  scatterer.  Ultrasonic  experiments  are  performed  for  two  different  geome¬ 
tries  of  three  submerged  rods.  To  aid  the  interpretation  of  the  results,  the  experiments  are 
numerically  modelled  by  the  multiple  multipoles  method. 

Introduction 

The  scattering  of  acoustic  waves  from  solids  submerged  in  a  fluid  has  received  well  deserved 
attention  for  a  long  time  (Everstine  and  Au-Yang,  1984;  Junger  and  Feit,  1986).  Even  for 
very  simple  cases,  the  scattered  wavefields  become  very  complex.  For  example,  in  the  case 
of  several  submerged  solid  objects,  the  major  events  can  readily  be  interpreted  as  direct 
reflections.  However,  the  remaining  wavefield  consists  of  comphcated  interactions  of  the 
wavefields  with  the  scatterers  and  the  water  surface.  Forward  modelling  of  the  situation  is 
a  valuable  aid  to  understand  these  interacting  events  better.  Especially  useful  are  methods 
which  decompose  the  scattered  wavefields  by  their  origin.  Thus,  forward  modelling  schemes 
such  as  the  finite  differences  method  have  to  be  discarded  as  interpretative  tools  because  they 
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only  yield  the  total  wavefields.  Various  other  methods  have  been  applied.  For  a  few  simple 
geometries,  analytic  solutions  are  available,  e.g.  cylindrical  objects  (Pao  and  Mow,  1973).  For 
more  complicated  geometries,  other  methods  have  been  used:  integral  equations  (de  Hoop, 
1990;  Luke  and  Martin,  1995),  the  T-matrix  (Bostrom,  1980),  perturbation  schemes  (Norris, 
1990),  or  boundary  element  methods  (Jensen  et  al,  1994)  which  might  be  coupled  to  finite 
elements  (Mathews,  1986;  Everstine  and  Henderson,  1990;  Fyfe  et  ai,  1991). 

Multiple  multipole  expansions  (MMP)  have  already  been  shown  to  be  a  versatile  tool  to 
solve  scattering  problems  in  fluids  (Imhof,  1995a)  or  solids  (Imhof,  1995b).  In  either  case, 
homogeneous  scatterers  were  embedded  in  a  homogeneous  fuUspace.  In  the  present  paper, 
a  combination  of  the  acoustic  with  the  elastic  scheme  is  introduced  as  an  interpretative  aid 
to  multiple  scattering  between  solids  submerged  in  a  fluid.  There  are  various  reasons  to 
use  MMP  expansions.  First,  they  have  been  shown  to  converge  faster  than  more  traditional 
methods  (Imhof,  1995a,b).  The  resulting  scheme  is  very  naturally  applied  to  situations 
of  multiple  scattering  objects.  Most  important,  the  scattered  fields  are  by  construction 
decomposed  by  scatterer.  Although  we  will  not  exploit  it,  the  wavefield  in  each  solid  scatterer 
is  also  decomposed  into  P-  and  S-waves. 

The  paper  is  structured  as  follows.  In  an  ultrasonic  watertank,  we  measure  the  scattering 
between  three  rods  of  Incite  and  the  water  surface.  Even  though  the  situation  is  rather  simple, 
we  find  the  scattered  field  to  contain  complex  interactions  between  the  rods  and  the  surface. 
To  aid  the  interpretation,  we  derive  the  acousto-elastic  MMP  scheme  which  we  use  to  model 
the  two  ultrasonic  experiments.  A  comparison  of  the  experimental  and  numerical  data  not 
only  validates  the  acousto-elastic  scheme  but  also  explains  various  events  as  internal  multiples 
in  the  rods  or  reverberations  between  the  scatterers  and  the  water  surface.  The  numerical 
data  also  identifies  some  events  as  reflections  being  diffracted  by  the  other  scatterers. 
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Ultrasonic  Modelling 


A  simple  fluid-solid  scattering  experiment  is  performed  in  an  ultrasonic  watertank.  Three 
rods  of  lucite  are  submerged  in  the  tank.  Two  different  geometries  are  used.  The  first  shown 
in  Figure  1  resembles  a  syncline  structure.  The  second  geometry  resembles  an  anticflne  shown 
in  Figure  2.  For  both  geometries,  an  acoustic  source  is  placed  above  the  rods.  Receivers 
placed  along  a  line  perpendicular  to  the  rods  measure  the  signals  reflected  from  the  rods. 
Although  the  ultrasonic  tank  is  of  finite  dimensions  (100  x  60  x  50  cm),  the  experiment  is 
designed  that  reflections  from  the  sides  and  the  bottom  of  the  tank  are  outside  the  time 
window  of  interest.  The  tank  is  equipped  with  a  PC-based  control  and  data  acquisition 
system.  A  schematic  thereof  is  presented  in  Figure  3.  Computer  controlled  holders  allow 
movement  of  a  source  and  a  receiver  along  a  prescribed  path  within  the  tank.  For  a  fixed 
source  -  scatterer  -  receiver  geometry,  the  computer  controlled  acquisition  system  allows 
to  perform  the  same  experiment  a  number  of  times  to  improve  the  signal-to-noise  ratio. 
Without  waveform  averaging,  the  smaller  events  would  vanish  in  the  ambient  noise  which 
reaches  the  same  amplitude  level  as  the  reflections. 

Piezoelectric  transducers  act  as  source  and  receiver.  The  source  (Panametrix  V323)  has 
a  strong  nonuniform  vertical  radiation  pattern,  as  most  of  its  energy  is  pointed  vertically 
downward  as  expected  from  a  vertical  point  force.  The  receiver  (Panametrix  V1091)  has  a 
receiver  pattern  very  similar  to  the  source  radiation  pattern.  This  source-receiver  combina¬ 
tion  tends  to  emphasize  waves  propagating  in  the  vertical  direction  and  to  suppress  waves 
propagating  in  the  horizontal  direction. 

The  rods  have  a  P-wave  velocity  a  =  2570  m/s,  a  S-wave  velocity  /?  =  1200  m/s,  and  a 
density  ps  of  1180  kg/m^.  The  P-wave  velocity  7  in  the  water  is  1460  m/s,  the  density  Pf 
is  1000  kg/m®.  The  width  of  the  first  scatterer  is  45mm,  the  second  60mm,  and  the  third 
30  mm.  All  scatterers  are  20  mm  thick.  Exiting  the  source  transducer  with  a  sharp  pulse,  the 
transducer  emits  a  wavelet  of  approximately  190  kHz  center  frequency  shown  in  Figure  4.  In 
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the  water,  the  resulting  dominant  wavelength  is  around  6  mm.  Thus,  the  scatters  are  3  to 
10  times  larger  than  the  propagating  pulse  which  allows  to  discriminate  between  reflections 
from  the  top  or  bottom  of  the  scatterers.  The  sampling  interval  is  chosen  to  be  400  ns  which 
yields  a  corresponding  Nyquist  frequency  of  1.25  MHz.  The  cutoff-frequencies  of  the  bandpass 
filter  are  set  to  5  KHz,  respectively  300  KHz.  For  each  trace,  512  samples  are  recorded.  The 
recording  is  delayed  by  120  //s  to  mute  the  direct  arrivals  and  to  maximize  the  time  window 
of  interest.  The  traces  are  averaged  over  2048  sweeps  to  reduce  the  noise  amplitudes.  The 
receiver  transducer  is  placed  at  75  different  positions  along  a  line  perpendicular  to  the  rods. 
The  spacing  between  the  receiver  points  is  5  mm.  The  source  transducer  is  located  in  the 
center  of  the  receiver  spread  at  position  38. 

Figure  5  shows  the  measured  record  for  the  synchne  model  defined  in  Figure  1.  As 
expected  for  a  syncline  geometry,  we  obtain  the  typical  butterfly  pattern  (Haartsen  et  al, 
1994).  Figure  6  shows  the  measured  record  for  the  anticline  model  defined  in  Figure  2.  For 
both  geometries,  the  major  events  can  readily  be  interpreted  as  direct  reflections.  But  the 
remaining  wavefields  consist  of  complicated  interactions  of  the  wavefields  with  the  scatterers 
and  the  water  surface.  Despite  the  simple  geometries,  it  is  not  obvious  how  to  interpret  the 
different  events.  As  a  consequence,  we  resort  to  forward-modelling  as  an  interpretative  aid. 

Theoretical  Background 

We  need  to  model  how  the  acoustic  pulse  described  by  the  displacement  u*"°(x,  u)  propagat¬ 
ing  in  the  water  scatters  from  submerged  solids.  To  distinguish  the  different  regions,  we  will 
use  the  symbol  F*^.  The  fluid  domain  is  denoted  by  F®.  Often,  the  fluid  region  will  also  be 
called  ‘background’.  The  three  solid  scatterers  are  denoted  by  F^,  F^,  and  F^.  The  boundary 
between  the  fluid  F®  and,  e.g.,  the  solid  F^  is  denoted  by  5Foi. 

In  the  frequency  domain,  the  displacement  u/  of  a  wave  travelling  in  a  two-dimensional. 
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homogeneous  fluid  is  described  by 


ivV-u/+u/  =  0,  (la) 

kf 

where  we  defined  the  wave  number  A;/  =  w/t  for  a  particular  frequency  u  and  the  propagation 
velocity  7.  The  fluid  velocity  can  also  be  written  as  7  =  y/^flpf  where  A/  is  the  Lame 
parameter  and  pf  the  density  of  the  fluid.  Note  that  we  suppress  a  common  time  factor  e  . 
Elastic  P-SV  waves  travelling  in  a  two  dimensional,  homogeneous  solid  are  described  by 

^VV-Us-ivx  Vxu,  +  u,  =  0,  (lb) 

where  we  defined  the  wave  numbers  kg  — (jO Jot  and  Ig  —  loJ (3  for  a  particular  frequency  lo  and 
the  propagation  velocities  a  =  -^/A^  +  2iXg/ Ps  and  (3  —  yfp>sfps-  The  parameter  A^,  and 
[ig  denote  the  density  and  the  Lame  parameters  of  the  medium. 

Instead  of  working  with  the  vector  forms  (1)  of  the  wave  equations,  we  will  use  the 
potentials  H(x,a;),  $(x,a;),  and  ^(x,u;)  =  y^(x,u;)  which  relate  to  the  displacements  u/ 


and  Uj  by 

u/(x,a;)  =  u“(x,  w)  =  VE(x,  w) ,  (2a) 

Us(x,6<;)  =  u^(x,a;) +  u^(x,a;)  =  V$(x,  w)  +  V  X '®'(x,a;) .  (2b) 

Each  scalar  potential  itself  satisfies  a  scalar  wave  equation: 

V^E(x,  a;)  +  A:^E(x,  a;)  =  0,  (3a) 

V^$(x,6c;)  +  A:^$(x,a;)  =  0,  (3b) 

V^^(x,6i;)  +  lf^(x,u;)  =  0.  (3c) 


Similar  to  the  purely  acoustic  or  elastic  cases  (Imhof,  1995a, b),  we  expand  the  potential 
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fields  due  to  the  D  —  3  scatterers  as 

D  Pd 

■^(x,  Cu)  =  ^  ^  ^  ^  ^  0,pnd  '  ^pndip^i  ^pdi  ^ f')  d”  ®  j 

d=l  P=1  n=-Wp^ 

$d(x,a;)  =  X!  XI  ^pn<i-<^pnd(x,x/d,A:,)  +e^,  (4b) 

P=1  n=-;v* 

^d(x,a;)  =  X  X  (h>nd-'tppnd{x,xldJs)  +el ,  (4c) 

P=1 

where  the  subscript  d  <  D  denotes  the  index  of  the  scatterer.  The  expansions  ^pn<i(x,  x^,  kf), 
</>p„d(x,x^,A:5),  and  ^pn<i(x,Xp^,^s)  are  all  solutions  to  their  respective  Helmholtz  equations 
(3).  For  each  scatterer  d  and  each  expansion,  e.g.  we  choose  a  set  of  P/  expansion  centers 
located  at  x^.  At  each  expansion  center  x^,  we  place  a  multipole  0pnd(x,  x^,  fc^)  of 
order  — A/pd  <  n  <  +Npd.  Because  the  maximal  orders  of  the  multipoles  are  finite  and  the 
resulting  expansion  is  non-orthogonal,  we  also  need  to  include  an  error  term  ef.  The  other 
two  potential  fields  E  and  ^  are  defined  in  exactly  the  same  way.  This  form  of  the  expansion 
allows  to  choose  different  numbers  and  locations  of  the  expansion  centers  for  each  scatterer 

and  each  field  E,  and  Also,  the  maximal  orders  for  the  different  multipoles  can  be 

chosen  independently. 

An  expansion  of  the  form  (4)  is  known  as  multiple  multipole  (MMP)  expansion  because 
we  have  in  fact  a  summation  over,  e.g.,  P/  multipoles  located  at  x^.  In  the  fluid,  we  choose 
the  propagatory  solutions  H^^{kfr)  to  the  Helmholtz  equation  (3a); 

^p„d(x,  x“d,  kf)  =  Hl^likf  |x  -  Xp^J)  with  iiiside  scatterer  d ,  (5a) 

where  is  the  angle  Z(x,  x  —  Xpd)  with  respect  to  the  unit  vector  x  in  the  x-direction. 
Foreseeing  the  need  of  a  free  surface  in  the  fluid  phase,  we  directly  add  the  stationary  phase 
contribution  of  the  field  reflected  at  the  surface  (Chew,  1988;  Imhof,  1995a)  to  the  expansion 
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function  (5a): 


,  k,)  =  H\^(k,  |x  -  xjl)  -  H\^(k,  |x  -  xfj)  e-'”-*  ,  {5b) 

where  =  Z(x,x  -  a;^).  The  location  of  the  corresponding  mirror  source  is  obtained 
by  reflection  of  the  expansion  center  x^  across  the  free  surface. 

In  the  solid,  we  can  choose  between  two  different  sets  of  expansion  functions  for  4>pnd  and 
ippnd-  Either  we  use  the  propagatory  solutions  and  or  we  choose 

the  standing  wave  solutions  J\n\{ksf’)  and  J|n|(^s’') 

if  x^  is  inside  scatterer  d  ,  (6a) 

if  x^  is  outside  scatterer  d  .  (6b) 


J\n\{ks  |x  -  x^rfl) 


4'pnd(p^i  ^pd^  ks) 


<‘>(fc.|x-x^|)e“V. 


in9^ 


‘H 


Similarly,  we  have: 


4)  = 


)ff<;>(4|x-x;j)e"-* 


if  x^  is  inside  scatterer  d  , 
if  x^  is  outside  scatterer  d  . 


(7a) 

(7b) 


The  three  different  wavefields  are  coupled  by  the  boundary  conditions  along  the  inter¬ 
face  of  the  fluid  and  the  solid.  The  boundary  conditions  require  continuity  of  the  normal 


displacement  n  •  u,  continuity  of  normal  traction  n  a  n  and  vanishing  tangential  traction 
t  •  (T  •  n.  For  a  point  x  on  the  interface  with  the  unit  normal  n  and  the  unit  tangential 
direction  t,  we  have 


n- (^u*”‘=(x,a;)-l-u-(x,a;))  =  n  •  (u^(x,a;)  H- u'*'(x,a;))  (8a) 

n  •  ^o-*”'^(x,a;)  -i-  <T“(x,a;))  •  n  =  n  •  ^<T*(x,a;)  +  a^{x,ujfj  •  n  (8b) 

0  =  t- ^cr*(x,a;)  +  <T'®'(x,a;)^ -n  (8c) 

where  the  displacements  are  defined  in  (2).  The  stress  tensor  in  the  fluid  is  given  by 
£r=(x,a;)  =  Ip(x,6c))  =  -Ifc^A/E:(x,a;)  which  also  defines  the  pressure  p(x,a;).  The  stress 
tensors  a^(x,  u>)  and  <t^(x,  w)  are  defined  in  a  local  coordinate  system  centered  at 
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the  expansion  center  of  the  multipole. 


-Xskl^  +  2/i 


Qj.2 


\  r  89“^  ) 


'1  1  8^' 
r  8r89  89  _ 


^4 


\r  89  ) 


^  -0 

[r  (,r  5^ 

>_o  J_^_ 

^  [2r2  89^ 


8^^V 
8r89)_ 
r  8  fl  8^  '\ 
2  8r  \r  5r  y 


(9) 

(10) 

(11) 


Note  that  the  components  of  displacements  and  of  the  stress  tensors  have  to  be  transformed 
or  rotated  into  the  global  coordinate  frame  (a:,  y,  z)  in  which  the  unit  normal  n  and  the  unit 
tangential  direction  t  are  defined. 

We  solve  for  the  unknown  coefficients  apnd,  I’pnd  and  Cpnd  by  enforcing  the  boundary 
conditions  (8)  on  discrete  matching  points  Xm  along  the  boundaries  of  the  scatterers.  We 
obtain  a  hnear  system  of  equations 


where  we  used  the  submatrices  S"^,  and  to  denote  the  normal  displacements  n  •  u 
at  the  matching  points  along  the  scatterer  s  G  {1,2,3}  due  to  ^pnd,  <i>pnd,  and  'ippnd-  The 
submatrices  S”^,  ,  and  are  the  same  but  for  the  normal  stress  n  •  <r  •  n.  The 

submatrices  and  contain  the  tangential  stresses  t  •  <r  •  n.  For  the  sake  of  clarity. 


the  index  d  is  used  to  indicate  which  scatterer  induces  the  field.  The  index  s  indicates  along 
which  scatterer  the  boundary  conditions  are  evaluated  and  thus  on  which  boundary  the 
matching  points  Xm  lie.  All  three  scatterers  contribute  to  the  scattered  field  H(xT„,a;)  in  the 
fluid.  Therefore,  we  have  submatrices  Ssd  for  each  scatterer  d  and  each  boundary  s.  But 
because,  e.g.,  scatterer  1  has  no  common  boundary  with  scatterer  2,  the  fields 
■^i{xm,oj)  and  ^2(xm,w)  do  not  interact  with  each  other.  Thus  for  s  #  d,  the 

submatrices  ^sd  and  ^sd  vanish.  The  individual  scatterers  are  only  coupled  by  the  scattered 
fields  Hi(x^,a;),  H2(x„„a;),  and  H3(x„^,a;)  propagating  in  the  fluid.  All  multiple  scattering 
is  automatically  accounted  for  by  the  boundary  conditions.  The  only  contribution  for  the 
fields  within,  e.g.,  scatterer  1  are  due  to  the  (self)  interactions  of  and  ^i(x,„,a;). 

Defining  the  matching  points  by  their  location  x^s  along  the  boundary  of  scatterer  s,  we 
can  write  the  submatrices  as 


^sd  —  ^^pn<i(xTns)  j  > 
^sd  ~  ^^pn<i(X77is:  j  ? 


=  [“n(0pnd(Xms,w))  , 

~  ^<^pnd(Xms?  ■> 

=  [^^nt(0pnd(Xms,a;))  , 


^sd —  ^'0pnd(Xms5  j  > 
^5d  =  [<^nn(^pnd(Xms,t^))]  , 
^"d=  [<^nt(^pnd(Xms,t^))]  • 

(13) 


Similarly,  the  vectors  and  o-””  on  the  right  hand  side  represent  the  incident  wavefield  in 
the  fluid  evaluated  at  the  matching  points  X7715. 


u?  = 

s 


(14) 


Finally,  the  vectors  ei<i,  and  cj,  contain  the  unknown  coefl3.cients  cipndi  ^pnd  and  Cpnd- 
Because  the  MMP  expansions  (4)  are  commonly  non-orthogonal,  the  system  (12)  is  made 
overdetermined  by  choosing  more  matching  points  than  needed.  It  is  then  solved  in  the 
least-squares  sense  by  QR  decomposition  (Strang,  1988)  using  Givens  updating  (Schwarz, 
1989). 
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Numerical  Modelling 


To  interpret  the  ultrasonic  results,  we  model  the  scattering  problems  depicted  in  Figures  1 
and  2.  Three  rods  of  Incite  are  submerged  in  water.  The  rods  have  a  P-wave  velocity 
a  —  2570  m/s,  a  S-wave  velocity  /5  =  1200  m/s,  and  a  density  Ps  of  1180  kg/m^.  In  the 
water,  the  P-wave  velocity  7  and  the  density  p/  are  1460  m/s,  respectively  1000  kg/m*.  A 
vertical  point  force  excites  an  acoustic  wave  which  propagates  downwards  and  interacts  with 
the  elastic  scatterers.  The  source  pulse  is  modulated  by  a  Ricker  wavelet  (Ricker,  1977) 
of  200  kHz.  The  resulting  scattered  fields  are  measured  at  75  locations  along  the  surface. 
Source,  scatterers  and  receivers  are  separated  by  at  least  30  dominant  wavelengths  of  6  mm. 
Due  to  the  resulting  long  propagation  times,  many  frequencies  are  needed  to  calculate  the 
traces.  For  a  sampling  interval  of  400  ns,  we  will  use  512  frequencies. 

To  damp  resonances,  we  account  for  the  intrinsic  attenuation  in  the  scatterers  by  making 
the  wave  number  k  complex 

k  =  kr  +  ia  (15) 


where  a  is  the  attenuation  coefficient.  Alternatively,  we  use  the  quality  factor  Q  which  relates 
to  the  attenuation  coefficient  a  by 


(16) 


where  the  velocity  c  is  either  a  or  j3.  Both  quality  factors  Qp  and  Qs  are  assumed  to  be  100. 

The  scattered  wavefield  in  the  fluid  is  expanded  as  (4a)  where  we  choose  £>  =  3,  P-^  =  4, 
P2"  =  5,  and  —  3.  For  each  of  these  multipoles,  we  choose  =  6.  As  expansion 
function,  we  use  (5b)  which  includes  the  reflections  from  the  free  surface.  Figure  7  deflnes 
the  exact  locations  of  the  expansion  centers  for  the  syncline  model  (Figure  1).  For  the 


anticline  model,  the  same  locations  relative  to  the  solids  are  used  for  the  expansion  centers. 

The  wavefields  in  the  solid  scatterers  are  expanded  as  (4b)  and  (4c)  for  the  P-,  respectively 
the  S-wave.  For  scatterer  we  choose  P^  =  P^  =  12.  For  scatterer  P^,  we  use  P2  = 
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=  14.  Finally,  for  the  third  scatterer  F®  we  employ  =  10  expansion  centers. 

At  each  expansion  center,  we  use  only  the  orders  between  —3  and  3.  Thus,  —  Np^  —  3. 
Because  we  placed  the  expansion  centers  for  the  wavefields  in  the  solids  outside  the  scatterers 
(Figure  7),  we  have  to  use  (6b)  and  (7b)  to  describe  the  P-,  respectively  the  S-waves.  For 
resonating  geometries  such  as  the  posed  problem,  the  propagatory  solutions  (6b)  and  (7b) 
are  superior  to  the  standing  wave  solutions  (6a)  and  (7a).  First,  their  amplitudes  decay 
faster  (Hafner,  1990).  But  they  also  force  us  to  illuminate  the  scatterers  from  all  sides  which 
decouples  the  multipoles  and  the  interfaces  because  each  multipole  illuminates  mainly  the 
interface  closest  to  its  expansion  center.  This  effect  lowers  the  condition  number  and  thus 
enhances  the  stability  of  the  numerical  matrix  inversion.  For  both  models,  660  expansion 
functions  and  444  matching  points  (=  1332  equations)  are  used. 

Figures  8  and  13  show  the  synthetic  seismograms  as  calculated  by  the  MMP  expansions 
for  the  S3mcline-,  respectively  the  anticline  model.  The  first  120  ^s  are  suppressed  since  they 
only  contain  the  direct  arrival. 

Comparison  of  the  Tank  Data  with  the  MMP  Solutions 

The  advantage  of  the  MMP  solution  over  the  ultrasonic  watertank  data  is  that,  by  con¬ 
struction,  the  wavefields  are  decomposed  by  scatterer.  Once  the  equation  system  (12)  is 
solved  for  a  particular  model,  we  fix  the  summation  index  d  to  1,  2,  or  3  in  expression  (4a) 
when  evaluating  the  seismograms.  Remember  that  the  index  d  denotes  the  scatterer  which 
emits  a  particular  field.  Thus,  we  can  plot  the  scattered  wavefields  for  each  of  the  three 
scatterers  independently  which  simplifies  the  seismograms  and  allows  correlation  of  events 
with  scatterers. 

In  the  ultrasonic  experiment,  we  use  a  point  source  and  a  point  receiver.  The  models 
contain  no  variation  perpendicular  to  the  source-receiver  plane.  Hence,  the  experiments  are 
in  fact  2|-D  experiments.  However,  the  numerical  MMP  scheme  is  derived  and  applied  to 
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2-D  only.  Comparing  the  experimental  results  to  the  numerical  ones  might  be  problematic. 
But  Esmersoy  (1986)  and  Lo  (1987)  showed  that  the  differences  between  2-D  and  2|-D  are 
neghgible  as  long  as  the  scatterers  are  in  the  far  field  with  respect  to  both  the  source  and 
the  receivers.  At  the  center  frequency  of  190  kHz,  the  scatterers  are  at  least  15  wavelengths 
away  from  either  source  or  receivers  for  both  geometries.  Clearly,  the  far  field  condition  is 
satisfied. 

Syncline  Model 

A  comparison  between  the  ultrasonic  record  (Figure  5)  and  the  MMP  solution  (Figure  8) 
shows  that  the  arrival  times  and  amplitudes  for  different  events  match  very  well.  The  main 
difference  between  the  Figures  5  and  8,  the  real  and  the  synthetic  data,  are  the  strong, 
repetitive  events  in  the  central  part  of  the  real  seismogram.  The  cause  is  the  signature  of 
the  source  tranducer  (Figure  4)  which  rings  more  than  a  Ricker  wavelet. 

The  butterfly  pattern  typical  for  synclines  can  easily  be  seen  (Haartsen  et  al,  1994).  As 
expected,  each  scatterer  reflects  twice  with  opposing  polarity  corresponding  to  reflections 
from  the  top  and  from  the  bottom.  The  two  events  are  typically  14  separated.  Toward 
the  end  of  the  traces  at  300  fis,  reverberations  of  the  reflection  from  scatterer  F^  can  be  found 
which  bounced  between  the  surface  and  the  scatters  F^  and  F^.  These  effects  can  be  seen 
more  clearly  in  the  MMP  solution  decomposed  by  scatterer  shown  in  Figures  9  to  11.  All 
synthetic  seismograms  are  scaled  similarly  which  allows  to  compare  the  amplitudes  between 
the  different  figures. 

The  decomposed  MMP  solution  in  Figures  9  to  11  shows  that  the  event  around  180  ns  in 
the  tank  data  (Figure  5)  is  in  reality  composed  of  different  events.  First,  there  is  the  direct 
reflection  from  scatterer  F^.  But  this  reflection  was  also  diffracted  through  the  scatterers  F^ 
and  F^  as  the  seismograms  9  and  11  show.  Also,  the  first  multiple  reflections  from  inside 
F^  and  F^  appear  nearly  at  that  time.  Another  interesting  set  of  event  appears  at  256  /iS 
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and  270  /iS.  Figure  9  correlates  these  two  events  with  scatterer  F^.  The  traveltimes  prove 
them  to  be  reflections  from  bouncing  between  the  water  surface  and  F  .  The  prior  one 
reflected  twice  from  the  top  of  F^  The  latter  one  reflected  once  from  the  top,  once  from 
the  bottom.  Thus,  one  would  expect  another  event  about  14  {is  later  corresponding  a  wave 
bouncing  between  the  bottom  of  F^  and  the  water  surface.  Indeed,  a  very  weak  event  appears 
at  the  expected  time  in  the  tank  data  (Figure  5)  and  the  MMP  solutions  (Figures  8  and  9). 
A  final  subtlety  are  the  multiples  arriving  around  300  pis  which  justify  the  stationary  phase 
reflections  into  the  MMP  expansion  functions  (5b).  The  scattered  fields  in  Figures  9  to  11 
show  very  nicely  that  these  events  are  a  combination  of  different  waves.  First,  there  are  waves 
reflected  from  F^,  then  bounced  of  the  water  surface  and  finally  rebounded  from  either  F^  or 
F^.  But  we  also  encounter  the  opposite:  waves  reflecting  from  F^  or  F^,  bouncing  of  the  water 
surface  and  finaUy  rebounding  from  F^.  To  demonstrate  the  effect  of  the  multiple  scattering, 
we  calculate  the  acousto-elastic  response  for  each  scatterer  separately  and  subtract  them 
from  the  complete  MMP  solution  (Figure  8).  The  residual  wavefield  is  shown  in  Figure  12. 
The  amplitudes  in  both  figures  are  scaled  by  the  same  amount  to  ease  a  direct  comparison  of 
single  scattering  versus  multiple  scattering.  Clearly,  the  two  scatterers  F^  and  F^  speed  up 
the  reflection  from  F^.  The  remaining  events  are  waves  bouncing  between  the  three  scatterers 
and  the  water  surface  as  well  as  internal  multiples  sped  up  passing  through  F^  and  F®. 

Anticline  Model 

Both  the  ultrasonic  record  (Figure  6)  and  the  MMP  solution  (Figure  13)  clearly  show  the 
reflections  from  the  top  and  the  bottom  of  each  scatterer.  Not  surprisingly,  scatterer  F 
shadows  F^  as  the  rather  abrupt  change  in  the  reflection-amplitude  around  170 /is  demon¬ 
strates.  To  identify  a  few  other  events,  we  decompose  the  MMP  solution  by  scatterer.  The 
decomposed  seismograms  are  shown  in  Figures  14  to  16. 

For  example.  Figure  15  shows  that  part  of  the  ringing  main  reflection  around  130  ^s  is 
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caused  by  multiple  scattering  inside  T^.  The  bottom  reflection  even  splits  into  2  events 
forming  a  butterfly  pattern  which  repeats  every  14  fis.  At  170  /zs,  the  main  reflection  from 
scatterer  appears  in  the  scattered  field  of  T^.  This  event  can  be  interpreted  as  reflection 
from  but  propagating  through  Similarly  at  190 /xs,  Figure  15  shows  the  reflection  of 

being  diffracted  through  F^.  After  both  events,  internal  multiples  from  and  F^  show 
up  in  the  fields  scattered  from  F^  and  F^,  respectively.  Diffractions  of  these  multiples  passing 
through  F^  are  again  present  in  Figure  15. 

To  demonstrate  the  effect  of  the  multiple  scattering,  we  calculate  the  response  of  each 
scatterer  alone  and  subtract  it  from  the  complete  MMP  solution  (Figure  13).  The  residual 
wavefield  is  shown  in  Figure  17.  The  amphtudes  in  both  figures  are  scaled  by  the  same 
amount  to  allow  a  direct  comparison  of  single  scattering  versus  multiple  scattering.  For  the 
first  170  iJS,  the  events  emanated  by  scatterer  F^  are  the  same  for  single  or  multiple  scattering. 
Clearly,  the  reflection  from  scatterer  F^  separates  into  two  events  in  the  residual  seismogram 
shown  in  Figure  17.  Scatterer  F^  not  only  shadows  F^  but  also  speeds  up  the  arrival  time  of 
the  reflection  passing  through  the  fast  solid.  The  same  holds  for  some  of  the  later  events  in 
Figure  17  corresponding  to  internal  multiples.  The  remaining  events  are  generated  by  waves 
bouncing  between  the  three  scatterers  and  the  water  surface.  For  example,  the  reflections 
around  300  /xs  are  caused  by  multiple  scattering  between  F^,  the  water  surface,  and  either  of 
F^  or  F^ 


Summary 

Already  for  simple  geometries,  the  wavefield  scattered  from  solids  submerged  in  water  can  be 
rather  complex.  Furthermore,  the  scattered  wavefield  becomes  even  more  complicated  if  the 
water  surface  interacts  with  the  scatters.  Especially  the  syncline  experiment,  where  three 
Incite  rods  were  submerged  in  an  ultrasonic  watertank,  demonstrates  this  point.  Although 
the  main  reflections  are  simple  to  identify,  the  remaining  events  are  very  hard  to  interpret. 
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To  aid  the  interpretation,  we  forward- modelled  the  situations  numerically  using  the  mul¬ 
tiple  multipole  method  (MMP).  For  each  submerged  solid,  the  scattered  wavefield  induced 
in  the  fluid  is  expanded  into  a  MMP  expansion  (Imhof,  1995a).  Similarly,  in  each  solid  two 
MMP  expansions  are  used  to  describe  the  P-  and  S-components  of  the  elastic  wavefields. 
The  advantage  of  the  MMP  scheme  is  that  the  scattered  wavefields  are,  by  construction, 
decomposed  by  scatterer  as  well  as  by  mode  of  propagation.  Hence,  the  wavefields  scattered 
from  each  solid  can  be  plotted  independently  which  facilitates  the  correlation  of  particular 
events  with  individual  scatterers.  Curiously,  this  allows  to  distinguish  direct  reflections  from 
reflections  diflFracted  by  another  scatterer  along  their  path  of  propagation. 

In  conclusion,  we  found  that  acousto-elastic  multiple  multipole  expansions  are  a  versatile 
tool  either  to  directly  solve  scattering  problems  or  to  aid  in  their  interpretation.  This  work 
should  find  uses  in  a  variety  of  fields,  notable  in  geophysics  and  oil  exploration,  as  well  as  in 
naval  underwater  applications. 
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Figure  1;  ‘Syncline’  model:  Three  oval  rods  of  Incite  are  submerged  in  water.  The  source- 
transducer  (star)  generates  a  pressure  wave  which  scatters  between  the  rods  and  the  water 
surface  until  it  is  picked  up  by  the  75  receivers  (triangles).  All  dimensions  are  in  millimeters. 
The  respective  centers  of  the  scatterers  F^,  and  F^  are  located  at  (—42.5,  —99),  (0,  —136) 
and  (65,  -98).  Both  the  source  and  receiver  38  are  positioned  at  (8, 0).  The  spacing  between 
the  receivers  is  5  mm. 
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Figure  2;  ‘Anticline’  model  where  the  star  denotes  the  source  and  the  triangles  symbolize 
the  receivers.  All  dimensions  are  in  millimeters.  The  respective  centers  of  the  oval  scatterers 
r\  r^,  and  r®  are  located  at  (-42.5,  -135),  (0,  -99)  and  (65,  -138).  Both  the  source  and 
receiver  38  are  positioned  at  (8, 0).  The  spacing  between  the  receivers  is  5  mm. 
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Network 


Figure  3:  Block  diagram  of  the  computer-based  ultrasonic  data-acquisition  system. 
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Figure  4:  Wavelets  after  propagating  30  wavelengths;  measured  in  ultrasonic  watertank 
(solid),  Ricker  wavelet  used  for  numerical  models  (dashed). 
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Figure  5:  The  seismogram  for  the  syncline  model  defined  in  Figure  1  as  measured 
ultrasonic  watertank. 
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Figure  7:  MMP  discretization  for  the  syncline  model.  Shown  are  the  locations  of  the  expan¬ 
sion  centers  for  the  multipoles.  The  multipoles  inside  the  scatterers  govern  the  wavefield  in 
the  fluid  background.  Multipoles  are  placed  all  around  the  scatterers  to  model  the  P-  and 
S- waves  on  the  inside.  The  discretization  for  the  anticline  model  is  the  same  but  upside-down. 
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Figure  9:  Syncline  model:  The  scattered  field  emanated  by  The  events  around  180 /iS 

are  reflections  originating  at  scatterer  but  passing  through  scatterer  F^  The  events  at 

300  fxs  are  reflections  from  F^  bouncing  between  the  surface  and  F^.  The  symbol ‘T’  denotes 

a  reflection  from  the  top  of  the  following  scatterer,  ‘B’  reflection  from  the  bottom,  ‘W’ 

reflection  from  the  water  surface,  and  ‘P’  meg.'pis  passes  through  the  scatterer. 
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Figure  10:  Syucline  model:  The  scattered  field  emanated  by  into  the  fluid. 
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Figure  11:  Syncline  model:  The  scattered  field  emanated  by  F^.  The  events  around  180 /iS 
are  reflections  from  scatterer  passing  through  scatterer  F®.  The  events  at  260 /iS  are 
reflections  from  F^  bouncing  between  the  water  surface  and  scatterer  F^.  The  events  at 
300  /xs  are  reflections  from  scatterer  F^  bouncing  between  the  surface  and  scatterer  F®. 
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Figure  15:  Anticline  model:  The  scattered  field  emanated  by  into  the  fluid.  The  events 
at  170  fxs  and  190  ^ls  are  the  reflections  from  F^  respectively  F^,  being  transmitted  through 
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Abstract 


We  analyze  the  influence  of  cavity  shape  on  far-field  seismic  wave  radiation  from  small 
explosion  sources  in  isotropic  and  transversely  isotropic  media  using  the  indirect  boundary 
element  method  (BEM).  The  analysis  utilizes  two  possible  mathematical  models  for  the 
pressure  fields  generated  by  the  explosion  in  the  acoustic  medium  inside  the  cavity.  One 
simple  model  for  the  source  is  a  volume  injection  (explosion)  point  source  inside  the  cavity, 
which  generates  a  pressure  field  that  arrives  approximately  instantaneously  at  each  point 
on  the  cavity  wall.  The  amplitude  of  the  field  decays  at  a  rate  inversely  proportional  to 
propagation  distance.  We  use  this  model  to  compare  results  with  finite  differences  solutions, 
and  the  test  confirms  the  accuracy  of  the  BEM  implementation.  Strong  compressional  and 
shear  signals  propagate  into  the  far-field  from  the  source  cavity,  and  the  compressional 
wave  amplitude  depends  strongly  on  direction  of  propagation.  The  behavior  is  similar  at 
frequencies  up  to  200  Hz,  including  those  typical  of  regional  seismic  wave  propagation  (1  to 
10  Hz). 

A  second  model  of  the  source,  one  that  has  been  applied  in  previous  analytical  and 
numerical  work,  is  to  assume  that  the  explosion  wiU  instantaneously  and  uniformlypxessmize 
the  source  cavity.  This  approach  yields  a  radiation  pattern  that  is  very  different  from  the 
point  source  model.  The  compressional  wave  radiation  patterns  have  a  large  variation  in 
signal  strength  and  very  strong  shear  waves  when  frequencies  up  to  200  Hz  are  included, 
but  seismograms  computed  for  frequency  ranges  of  interest  in  regional  wave  propagation 
(around  1  Hz)  have  a  nearly  isotropic  compressional  wave  radiation  pattern  and  a  much 
smaller  shear  wave.  This  shear  wave  is  largest  for  long,  tunnel-like  cavities  or  short,  disk-like 
cavities.  When  the  cavity  is  located  in  a  transversely  isotropic  medium,  quasi-shear  waves 
are  generated,  but,  for  the  media  we  consider,  their  magnitude  decreases  for  long,  tunnel-like 
cavities  and  is  larger  for  more  equidimensional  cavities.  These  somewhat  counterintuitive 
results  show  that  anisotropy  can  actually  cause  a  radiation  pattern  to  appear  more  isotropic 
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than  the  corresponding  wavefields  in  an  isotropic  medium. 


Introduction 

Analysis  of  seismic  waves  generated  by  nuclear  explosions  is  complicated  by  several  factors. 
For  example,  spall  and  tectonic  release  processes  contribute  strongly  to  far-field  seismic  wave 
signals,  and  the  classic  explosion  source  model,  which  yields  an  isotropic  compressional  wave 
radiation  pattern,  cannot  reproduce  observed  waveforms  (e.g.,  Wallace  et  al.,  1983,  1985; 
Priestley  et  al,  1990;  Patton,  1991;  Schlittenhardt,  1991).  Changes  in  source  cavity  shape 
and  source  location  within  the  cavity  also  can  influence  the  radiation  pattern  of  the  explosion. 
Zhao  and  Harkrider  (1992)  developed  solutions  for  the  displacement  field  created  by  an  ex¬ 
plosion  source  located  off-center  within  a  soUd  sphere  embedded  in  an  infinite  solid  medium, 
a  model  of  a  fully  tamped  explosion.  Here  the  loss  of  symmetry  causes  shear  waves  to  be 
generated,  even  with  the  spherical  cavity.  Other  analyses  have  examined  instead  the  effect  of 
varying  the  cavity  shape.  Glenn  et  al.  (1985,  1986),  Rial  and  Moran  (1986)  and  Glenn  and 
Rial  (1987)  considered  ellipsoidal  and  cylindrical  cavity  shapes,  and,  applying  a  combination 
of  numerical  and  approximate  analytical  methods,  obtained  estimates  of  far-field  radiation 
patterns.  They  suggested  that  the  P-wave  radiation  would  be  strongest  perpendicular  to  an 
ellipsoidal  cavity  and  that  a  fairly  strong  S-wave  would  be  generated  as  the  cavity  aspect 
ratio  (ratio  of  length  to  diameter)  became  large.  Stevens  et  al.  (1991)  also  examined  an 
ellipsoidal  cavity.  However,  in  contrast  with  the  earher  results,  they  concluded  that  at  the 
frequencies  of  interest  for  regional  wave  propagation,  the  P-wave  radiation  pattern  would  be 
fairly  isotropic  and  S-wave  amplitudes  comparatively  small. 

We  continue  these  investigations  by  using  the  indirect  boundary  element  method  (BEM) 
to  simulate  the  seismic  wavefields  generated  by  explosions  in  cylindrical  cavities.  The  BEM 
is  a  method  for  the  computation  of  full  waveform  synthetic  seismograms  that  discretizes 
boundaries  in  an  earth  model  into  elements  and  then  replaces  them  with  fictitious,  secondary 
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sources  (Chew,  1990).  It  has  recently  been  applied  to  seismic  exploration  problems  where  a 
source  is  located  inside  infinite  and  semi-infinite  cylindrical  boreholes  (e.g.,  Bouchon,  1993; 
Dong,  1993;  Dong  et  ai,  1995;  Dong  and  Toksoz,  1995).  One  major  advantage  of  this 
approach  is  that  unlike  finite  difference  schemes  it  does  not  require  discretization  of  the 
entire  earth  model  and  can  therefore  be  used  to  compute  wavefields  at  large  distances  from 
a  comparatively  small  borehole.  We  have  now  extended  the  method  to  apply  to  a  cylindrical 
explosion  source  cavity  of  finite  length.  This  cavity  model  could  represent,  for  example,  a 
horizontal  tunnel  or  a  vertical  shaft  containing  an  explosive  device.  The  cavity  is  assumed  to 
be  filled  with  an  acoustic  medium,  and  it  is  embedded  within  an  infinite,  homogeneous  elastic 
medium.  This  elastic  material  can  be  either  isotropic  or  transversely  isotropic  with  the  axis 
of  symmetry  parallel  to  the  cavity  axis.  Once  velocities  and  elastic  constants  are  specified 
for  the  two  media,  the  BEM  gives  full  waveform  synthetic  seismograms  corresponding  to 
the  seismic  wavefields  produced  by  a  source  inside  the  cavity  and  receivers  in  the  elastic 
formation.  Assuming  the  the  explosion  is  a  small,  fully  decoupled  event,  we  can  apply  this 
linear  simulation  to  the  computation  of  far-field  synthetic  seismograms. 

After  describing  the  BEM  in  more  detail  below,  we  present  computational  results  using 
two  different  mathematical  models  for  the  cavity  source.  First,  we  examine  the  far-field 
seismic  waves  generated  by  a  volume  injection  point  source  located  on  the  axis  of  the  cavity. 
A  major  incentive  for  using  this  source  model  is  that  it  was  available  in  existing  finite  differ¬ 
ence  codes  and  therefore  allowed  a  straightforward  calibration  of  the  BEM.  A  comparison 
of  results  from  the  two  methods  shows  that  our  implementation  of  the  BEM  is  accurate. 
We  next  present  radiation  patterns  varying  both  cavity  shape  and  location  of  the  source 
within  the  cavity,  noting  that  the  strength  of  the  incident  pressure  field  depends  on  distance 
from  the  source.  Radiation  patterns  show  strong  variations  in  P-wave  amplitude  and  strong 
S-waves  at  all  frequencies  considered. 

The  second  source  model  we  consider  makes  the  assumption  that  the  explosion  instan¬ 
taneously  and  uniformly  pressurizes  the  cavity,  so  that  the  source  location  within  the  cavity 
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is  no  longer  a  relevant  parameter.  This  model  is  probably  a  more  realistic  description  of  the 
explosion  process,  and  it  has  been  applied  frequently  in  the  past  (e.g.,  Stevens  et  al.,  1991; 
Patton,  1991).  It  also  yields  radiation  patterns  that  are  quite  different  from  the  point  source 
model.  The  results  we  present  below  show  that  at  low  frequencies  (1  Hz)  the  source  creates 
a  more  nearly  isotropic  P-wave  radiation  pattern  and  that  the  S-wave  is  relatively  small, 
though  it  depends  on  cavity  aspect  ratio.  The  solutions  from  this  low  frequency  range  up 
to  200  Hz  are  similar  to  the  results  of  Stevens  et  a/.(1991),  and  they  suggest  that  placing 
a  source  inside  a  tunnel  instead  of  a  more  equidimensional  source  cavity  may  not  generate 
very  strong  variations  in  far-field  seismic  waves.  Continuing  with  the  instantaneous  pres¬ 
surization  model,  we  apply  it  to  the  computation  of  radiation  patterns  for  both  hard  and 
soft  rock  formations  to  show  the  dependency  of  the  wavefields  on  formation  properties.  In 
addition,  we  also  present  the  radiation  patterns  for  hard  and  soft  transversely  isotropic  me¬ 
dia.  The  solutions  for  these  media  demonstrate  that  the  quasi-shear  wave  generated  by  an 
explosion  source  can  actually  be  smaller  than  the  analogous  results  for  isotropic  formations 
when  the  cavity  has  a  large  aspect  ratio.  Hence  we  make  the  somewhat  surprising  prediction 
that  an  explosion  in  an  anisotropic  medium  may  generate  less  shear  energy  than  the  equiv¬ 
alent  source  located  in  an  isotropic  formation,  a  result  that  can  be  explained  by  considering 
equivalent  moment  tensor  representations  of  the  sources. 

Boundary  Element  Method 

The  indirect  boundary  element  method  (BEM)  has  recently  been  applied  to  the  calculation 
of  acoustic  and  elastic  wavefields  in  and  around  infinite  boreholes  in  both  homogeneous  and 
plane  layered  media  (Bouchon,  1993;  Dong,  1993;  Dong  et  al.,  1995;  Dong  and  Toksoz,  1995). 
In  addition,  the  method  has  been  extended  for  application  to  semi-infinite  cavities  to  model 
the  radiation  from  seismic  sources  on  or  near  the  bottom  end  of  a  borehole  (Dong  et  al, 
1995).  The  basic  idea  of  this  approach  is  that  the  effects  of  the  diffracting  boundary  on  the 
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displacement  fields  in  the  elastic  medium  can  be  represented  by  a  set  of  fictitious  sources 
located  on  the  boundary,  the  cavity  wall.  Once  these  fictitious  sources  and  the  Green’s 
functions  for  the  elastic  medium  are  known,  the  wavefield  generated  by  the  source  inside  the 
cavity  can  be  calculated  for  a  point  outside  the  cavity.  Our  work  extends  these  earlier  results 
to  include  the  second  end  of  the  cylindrical  cavity,  allowing  synthesis  of  the  displacement 
fields  generated  by  an  explosion  source  located  inside  a  tunnel. 

Considering  a  point  source  in  the  acoustic  medium  on  the  axis  of  the  cavity,  the  displace¬ 
ment  field  can  be  related  to  the  following  potential  (Dong  et  al,  1995): 

<^(x)  =  <{)i+  f  dS'g{yi,J(.')ip{x')  (1) 

JB 

Here  is  the  potential  created  by  the  point  source,  p(x,  x')  is  the  scalar  Green’s  function 
for  the  acoustic  medium,  and  '^(x^)  is  an  acoustic  source  distribution  along  the  cavity  wall. 
The  integrand  is  evaluated  on  the  cavity  surface  B.  Equation  (1)  yields  the  displacement 
field  u(x,  w)  or  pressure  field  p(x,  co)  anywhere  in  the  cavity  or  on  the  cavity  wall  via 

u(x,  uj)  =  V0(x,  w),  p(x,  a;)  =  -AV^0(x,  lo)  (2) 


Similarly,  the  displacement  field  U(x)  in  the  elastic  medium  surrounding  the  cavity  is  given 
by 


U(x)  =  dS"G(x,x')  •  ’$^(x'). 


(3) 


Here  G(x,x')  is  the  Green’s  tensor  for  the  medium,  which  in  our  implementation  can  be 
either  isotropic  or  transversely  isotropic  with  the  axis  of  symmetry  parallel  to  the  cylinder 
axis  (Dong  et  al,  1995).  The  vector  force  source  on  the  cavity  wall  is  denoted  by  '^(x') 
and  consists  of  two  components  in  the  vertical  and  radial  directions.  Because  of  the  axial 
symmetry  inherent  in  the  prescribed  geometry  of  the  source  and  cylindrical  cavity,  the  third 
component  of  the  force  vector  is  identically  zero. 

The  numerical  algorithm  begins  by  discretizing  the  wall  of  the  cylindrical  cavity  into  rings 
on  the  cylinder  and  into  circular  annuli  on  the  top  and  bottom  capping  surfaces  (Figure  1). 
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In  general,  there  must  be  at  least  3  elements  within  the  smallest  wavelength  involved  in  the 
calculation.  The  sources  are  assumed  to  be  constant  on  each  element,  and  the  magnitudes 
of  these  acoustic  and  force  sources  are  determined  by  enforcing  boundary  conditions  on 
each  element.  These  conditions  are  the  continuity  of  the  normal  components  of  stress  and 
displacement  and  the  vanishing  of  the  tangential  stress  in  the  solid  at  the  boundary  between 
the  acoustic  and  elastic  media.  By  writing  these  equations  in  matrix  form,  we  obtain  a 
system  which  can  be  solved  for  the  source  strengths  on  each  of  the  Ne  elements: 


Ne  Nc 


t=l  i=l  i=l 

W.  Ne  Ne 

(4) 

EB/^+E^Ji^i’  +  EBli^:''  =  O’" 

i=l  i=l  i=l 

Ne  Ne 

(5) 

j:q,Fr+Y:qiN  =  of- 

(6) 

Z=1  Z=1 


and  are  the  displacement  fields  at  the  yth  element  created  by  the  volume  in¬ 
jection,  vertical  force  and  radial  force  sources,  respectively,  on  the  zth  element.  They  are 
determined  by  integrating  the  Green’s  functions  over  each  element.  In  the  same  way,  the 
quantities  B  and  C  are  the  radial  and  tangential  stresses  on  each  element.  The  variables  D, 
following  similar  notation  conventions,  are  the  exciting  fields,  or  incident  wavefields  gener¬ 
ated  by  the  acoustic  source,  on  each  element.  There  are  a  total  of  3Ne  equations,  which  can 
be  solved  for  the  acoustic  sources  V/  and  force  sources  F[,  Fi  once  the  different  acoustic 
and  elastic  displacement  and  stress  fields  are  known.  These  are  computed  using  discrete 
wavenumber  methods  to  evaluate  integrals  over  horizontal  wavenumber  for  the  three  types 
of  media  involved  in  our  simulations  (acoustic,  isotropic  and  transversely  isotropic)  (Dong 
et  al,  1995). 

In  our  calculations,  we  have  utilized  two  different  models  of  the  exciting  wavefields  D. 
The  first  approach  was  to  consider  a  point  source  inside  the  cavity.  The  incident  pressure  and 
displacement  can  then  be  computed  using  a  discrete  wavenumber  integral  solution.  When 
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the  acoustic  velocity  in  the  cavity  is  very  high,  the  arrival  time  differences  of  the  wavefield 
along  the  cavity  wall  will  be  negligible,  though  the  amplitude  will  inversely  proportional  to 
distance  from  the  source.  A  second  model  of  the  source  field  is  to  assume  instantaneous 
pressurization  of  the  cavity  by  the  explosion  so  both  the  amplitude  and  the  arrival  time 
of  the  pressure  signal  are  constant.  This  assumption  is  commonly  made  in  the  analysis  of 
nuclear  explosions  in  spherical  cavities  (e.g.,  Patton,  1991).  To  apply  this  source  model,  we 
simply  set  DJ  =  =  0  and  =  1.0  as  the  source  condition.  Below  we  show  results 

for  both  source  models,  demonstrating  the  accuracy  of  the  BEM  and  also  highlighting  the 
difference  in  radiation  patterns  predicted  by  the  two  different  methods. 

Radiation  Pattern  Results 


Point  Source 

We  applied  the  indirect  BEM  to  point  sources  in  cylindrical  cavities  in  a  hard  rock  medium 
with  P-wave  velocity  4000  m/s,  S-wave  velocity  2200  m/s,  and  density  2200  kg/m^.  The 
acoustic  velocity  and  density  of  the  cavity  filling  material  were  set  to  10000  m/s  and 
0.129  kg/m^,  respectively,  to  model  a  rapid  pressure  wave  generated  by  an  explosion.  A 
comparison  of  BEM  results  with  finite  difference  simulations  provides  a  useful  demonstra¬ 
tion  of  the  accuracy  and  rapid  computation  speed  of  the  BEM,  and  we  implemented  this  test 
with  a  cylindrical  cavity  model  of  length  160  m  and  radius  40  m  (with  aspect  ratio,  the  ratio 
of  length  to  diameter,  therefore  equal  to  2).  Receivers  were  located  at  a  constant  distance 
of  1250  m  from  the  cavity  center,  these  locations  covering  the  entire  range  of  angles  from 
along  the  cavity  axis  to  perpendicular  to  this  axis  in  order  to  estimate  radiation  patterns. 
This  comparatively  small  distance  of  the  receivers  was  used  to  make  computation  time  with 
the  finite  difference  algorithm  practical,  since  a  relatively  fine  discretization  of  the  model 
is  required  to  accurately  model  the  source  cavity.  A  major  disadvantage  of  finite  differ- 
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ence  methods  for  modeling  source  behavior  is  this  discretization  requirement,  as  increasing 
source/receiver  distances  to  more  realistic  values  greatly  increases  computation  times.  Be¬ 
cause  this  distance  was  so  small,  we  also  set  the  center  frequency  of  the  source  time  wavelet 
to  8  Hz  to  obtain  far-field  radiation  patterns.  Comparisons  of  the  synthetics  seismograms 
from  the  two  algorithms  show  that  the  results  are  essentially  identical  (Figure  2),  but  with 
these  parameters,  the  finite  difference  simulation  required  about  13  hours  to  computed  1  sec 
time  series  on  a  Digital  Alpha  workstation,  compared  to  about  2  hours  45  minutes  for  the 
BEM.  The  latter  computation  computed  frequency  domain  results  up  to  34  Hz,  yielding 
synthetic  seismograms  3.75  sec  in  length.  For  this  modeling  task,  the  BEM  clearly  is  a  more 
useful  tool  than  finite  differences. 

Another  important  advantage  of  BEM  not  revealed  by  this  comparison  is  that  it  provides 
the  complete  solution  with  no  auxiliary  computations.  Other  approaches  to  modeling  radi¬ 
ation  from  source  cavities  couple  two  methods,  one  for  near  source  propagation,  another  for 
propagation  to  the  far-field.  For  example,  the  representation  theorem  can  be  used  to  esti¬ 
mate  far-field  solutions  after  an  initial  computation  with  finite  differences  or  finite  elements 
(e.g.,  Glenn  et  al,  1985;  Stevens  et  ai,  1991).  The  BEM  avoids  a  second  step  by  directly 
computing  the  wavefields  at  an  arbitrary  distance. 

Radiation  patterns  were  estimated  for  the  BEM  and  finite  difference  seismograms  by 
rotating  the  vertical  and  horizontal  traces  into  components  parallel  and  perpendicular  to  the 
propagation  direction  to  isolate  the  compressional  and  shear  waves,  respectively,  on  different 
traces.  The  amplitudes  were  then  estimated  by  taking  the  maxima  of  the  envelopes  of  the 
individual  traces.  While  there  are  some  slight  differences  in  these  radiation  patterns,  the 
results  are  clearly  very  close  and  this  test  confirms  that  the  BEM  works  accurately  (Figure  3). 

In  order  to  explore  the  dependence  of  the  radiation  pattern  on  cavity  size,  we  computed 
radiation  patterns  for  a  several  cavity  shapes  using  the  BEM  (Figure  4).  Each  cavity  model 
had  a  length  of  80  m,  and  the  radius  of  the  cavity  was  increased  in  order  to  decrease  the 
aspect  ratio.  The  source  was  at  the  center  of  the  cavity.  Receivers  were  located  at  a  distance 
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of  10000  m  to  examine  far-field  radiation  patterns,  and  the  frequency  domain  response 
was  convolved  with  a  wavelet  having  a  dominant  frequency  of  1  Hz.  As  the  aspect  ratio 
diminishes,  the  P-wave  pattern  tends  to  become  more  isotropic,  while  the  S-wave  amplitude 
decreases  relative  to  the  P-wave.  For  larger  aspect  ratios,  which  might  be  more  realistic 
models  of  tunnel  source  cavities,  the  S-wave  is  as  large  as  the  compressional  wave  signal. 

Variations  in  aspect  ratio  clearly  yield  significant  changes  in  radiation  pattern,  but  an¬ 
other  important  variable  for  the  point  source  explosion  model  is  the  location  of  the  volume 
injection  source  within  the  cavity.  Though  the  method  as  currently  implemented  requires 
the  source  to  be  on  the  axis  of  the  cavity,  we  are  free  to  move  it  anywhere  along  this  axis. 
In  Figure  5,  we  show  radiation  patterns  for  the  cavity  with  an  aspect  ratio  of  5  where  the 
source  is  located  at  several  different  points  between  the  center  of  the  cavity  and  one  end. 
As  the  source  approaches  the  end  of  the  cavity,  the  radiation  pattern  becomes  essentially 
identical  to  that  of  a  single  point  force  in  an  infinite  elastic  medium.  This  variation  in  the 
far-field  signals  generated  by  the  source  is  caused  by  the  variation  of  incident  amplitude 
with  distance  from  the  source.  Since  the  amplitude  is  inversely  proportional  to  distance,  the 
incident  field  at  the  far  end  of  the  cavity  is  almost  two  orders  of  magnitude  smaller  than 
that  at  the  near  end  when  the  source  is  located  1  m  from  the  end.  Hence,  the  reradiation  of 
energy  by  the  closer  end  of  the  cavity  dominates  the  wavefields. 

Instantaneous  Pressurization  Source 

Like  the  point  source,  the  instantaneous  cavity  pressurization  source  model  will  be  a  func¬ 
tion  of  cavity  aspect  ratio,  though  the  position  of  the  explosive  device  within  the  cavity  is 
obviously  irrelevant.  In  Figure  6  we  compare  the  spectra  from  0.1  to  200  Hz  of  the  far-field 
displacements  computed  with  the  volume  injection  point  source  and  cavity  pressurization 
models.  Spectra  were  computed  for  angles  of  3,  45  and  87  degrees  from  the  axis  of  symmetry 
in  order  to  avoid  slight  numerical  instabilities  inherent  in  the  method  for  the  angles  of  0  and 
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90  degrees.  The  results  in  Figure  6A  show  that  the  spectrum  for  the  pressurization  model 
displays  a  flat  response  out  to  a  corner  frequency  that  depends  on  the  direction  of  observa¬ 
tion.  This  corner  frequency,  around  50  Hz,  is  larger  for  propagation  perpendicular  to  the 
cavity  axis  than  along  that  axis.  In  addition,  most  of  the  difference  in  amplitudes  is  found 
at  frequencies  well  over  10  Hz.  Similar  results  were  obtained  by  Glenn  et  al.  (1985,  1986) 
and  Stevens  et  al.  (1991).  However,  the  point  source  spectrum  shows  a  linear  behavior  for 
frequencies  from  0.1  to  approximately  50  Hz,  (Figure  6B)  and  the  difference  in  amplitudes 
is  fairly  constant  throughout  the  frequency  range  displayed  in  the  figure. 

One  very  important  feature  of  the  instantaneous  pressurization  spectra  (Figure  6A)  is 
that  most  of  the  directional  variation  occurs  for  frequencies  larger  than  50  Hz.  The  radiation 
pattern  for  the  step  function  response  does  predict  a  large  variation  in  amplitude  of  the 
P-wave  with  observation  direction  relative  to  the  cavity  axis,  much  like  the  point  source 
(Figure  7).  The  shear  wave  amplitude  is  also  relatively  large  compared  to  the  compressional 
signal.  However,  when  the  synthetic  seismograms  are  produced  by  convolving  the  spectral 
domain  response  with  a  Ricker  wavelet  of  1  Hz  to  examine  the  lower  frequencies  that  are 
more  typical  of  regional  or  teleseismic  wave  propagation,  the  S-wave  strength  is  reduced 
compared  to  the  P-wave  and  the  compressional  wave  amplitude  is  nearly  independent  of 
propagation  direction  (Figure  8).  In  this  situation,  the  point  source  still  shows  a  strong 
shear  wave,  but  the  pressurization  source  model  now  displays  a  uniform  P-wave  radiation 
pattern  and  a  much  smaller  S-wave. 

Behavior  of  the  pressurization  source  is  further  demonstrated  in  Figure  9,  where  we 
show  compressional  and  shear  wave  radiation  patterns  obtained  by  convolving  the  frequency 
domain  displacement  field  with  a  Ricker  wavelet  with  a  1  Hz  center  frequency  for  several 
aspect  ratios.  Computations  for  this  figure  used  a  constant  cavity  radius  of  4  m  and  varied 
the  length  to  assess  the  effect  of  changing  aspect  ratio.  The  shear  wave  is  smallest  relative 
to  the  P-wave  for  aspect  ratio  near  1. 
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Low  Velocity  Formation 

The  radiation  pattern  of  the  pressurization  source  will  have  some  dependence  on  the  velocity 
of  the  surrounding  elastic  medium  as  well  as  on  the  shape  of  the  cavity.  Radiation  patterns 
for  a  soft  rock  formation  are  similar  to  those  obtained  for  the  hard  rock  formation  (Figure  10). 
The  P  and  S-wave  velocities  of  the  soft  rock  medium  are  2290  m/s  and  980  m/s,  respectively, 
and  the  density  is  2250  kg/m^.  These  velocities  are  approximately  those  of  the  Pierre  Shale, 
which  is  actually  transversely  isotropic  (White  et  ai,  1983).  The  principal  difference  between 
these  patterns  and  those  in  Figure  9  is  that  the  radiated  fields  in  the  shale  are  much  larger. 
If  we  were  to  represent  the  cavity  source  by  an  equivalent  moment  tensor  source,  the  far-field 
amplitude  of  the  displacement  field  would  be  inversely  proportional  to  the  cubed  velocity 
of  the  formation,  which  explains  the  larger  amplitudes  in  the  slow  medium  (Ben-Menahem 
et  al,  1991).  At  the  same  time,  the  interaction  of  the  cavity  with  the  medium  will  change, 
leading  to  a  different  equivalent  moment  tensor  solution.  The  largest  P-wave  amplitudes 
compared  to  S- waves  are  found  for  aspect  ratio  1.25  in  the  hard  rock  formation  and  for 
aspect  ratio  2.5  in  the  soft  rock.  Hence,  it  is  clear  that  the  equivalent  moment  tensors  in  the 
two  media  are  somewhat  different,  leading  to  subtle  differences  in  seismic  wave  radiation. 

Transversely  Isotropic  Media 


Aligned  Fractures 

A  medium  containing  aligned  cracks  is  a  realistic  example  of  a  material  that  is  effectively 
anisotropic  for  long  wavelengths  (Hudson,  1980,  1981).  This  transversely  isotropic  medium 
can  be  incorporated  in  the  BEM  calculations  when  the  axis  of  symmetry  is  parallel  to  the 
cylindrical  cavity  axis,  which  occurs  when  the  cracks  are  perpendicular  to  the  axis.  If  the 
cavity  represents  a  horizontal  tunnel,  then  the  aligned  fractures  would  all  be  vertical. 

We  estimated  the  elastic  constants  of  such  a  material  using  the  theory  of  Hudson  (1980, 
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1981),  in  which  the  effective  elastic  constants  of  the  anisotropic  medium  are  given  in  terms  of 
the  isotropic  parameters  of  the  unfractured  rock  plus  perturbation  terms  accounting  for  the 
effects  of  the  aligned  cracks.  We  used  the  velocities  of  the  isotropic  material  from  the  radia¬ 
tion  pattern  computations  shown  above  for  the  background  parameters.  The  perturbations 
were  computed  using  a  crack  density  ^  =  na^  =  0.05  (n=number  density  of  cracks,  a=radius 
of  the  penny-shaped  cracks)  and  assuming  dry  (empty)  cracks.  Table  1  displays  the  elastic 
constants  of  both  the  isotropic  material  and  the  composite  medium  with  cracks,  while  Fig¬ 
ure  11  shows  the  phase  velocity  variations  of  the  three  wave  types,  the  quasi-compressional 
wave  (qP),  the  fast  quasi-shear  wave  (qSa)  and  the  slow  quasi-shear  wave  (qSi).  If  the  cavity 
axis  is  defined  as  the  2;-axis,  then  the  qSi-wave  is  polarized  in  the  vertical  plane  containing 
the  ray  from  source  to  observation  point  and  will  be  called  the  qSV  wave  for  simplicity. 

The  qP-wave  radiation  patterns  for  the  explosion  source  are  very  similar  in  magnitude  and 
shape  to  those  obtained  for  the  isotropic  formation  (Figure  12).  However,  the  qSV  radiation 
patterns  show  different  amplitude  variations  than  the  isotropic  case.  Figure  13  emphasizes 
this  point  by  directly  comparing  radiation  patterns  for  the  isotropic  and  anisotropic  media 
at  aspect  ratios  of  5  and  1.25.  The  most  interesting  point  is  that  the  minimum  qSV  am¬ 
plitude  is  found  for  the  largest  aspect  ratio  value  of  5,  and  the  relative  amplitude  increases 
monotonically  with  decreasing  aspect  ratio. 

This  result  is  at  first  somewhat  counter-intuitive,  since  computations  of  radiation  patterns 
using  the  classic,  point  source  explosion  model  (three  equal,  mutually  perpendicular  dipoles) 
show  comparatively  large  quasi-shear  wave  amplitudes  (Mandal  and  Toksoz,  1991;  Ben- 
Menahem  et  al,  1991).  However,  this  difference  can  be  easily  explained.  The  radiation 
patterns  for  the  isotropic  medium  (Figure  9)  correspond  to  the  far-field  radiation  of  a  source 
consisting  of  two,  equal  and  perpendicular  horizontal  dipoles  and  a  smaller  vertical  dipole 
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(e.g.,  Glenn  et  al,  1986;  Rial  and  Moran,  1986).  This  moment  tensor  has  the  form 

Dh  0  0 

M  =  0  D;,  0  ,  (7) 

0  0 

where  Dh  and  Dy  are  the  magnitudes  of  the  horizontal  and  vertical  dipoles,  respectively.  The 
relative  magnitudes  of  these  two  dipoles  for  the  isotropic  medium  are  straightforwardly  esti¬ 
mated  from  the  amplitudes  of  P-waves  propagating  in  the  horizontal  and  vertical  directions, 
since  all  other  factors  affecting  amplitudes  are  constant  with  direction.  The  amplitudes  of 
the  dipoles  can  also  be  measured  for  the  anisotropic  medium,  but  the  variation  in  geomet¬ 
rical  spreading  with  direction  must  be  considered  (Ben-Menahem  et  al,  1991).  Simple  ray 
tracing  procedures  can  be  used  to  compute  the  spreading  values  for  horizontal  and  vertical 
directions  (e.g.,  Gibson  et  oZ.,  1991),  showing  that  the  moment  tensors  are  almost  identical 
for  the  isotropic  and  fractured  media.  For  example,  when  cavity  aspect  ratio  is  5,  we  find 
Dy  =  0.9Dh.  Computation  of  ray  theoretical  radiation  patterns  with  this  moment  tensor  rep¬ 
resentation  confirm  that  the  some  moment  tensor  yields  the  radiation  patterns  in  Figures  9 
and  12. 

This  moment  tensor  representation  of  the  source  allows  a  simple  explanation  of  the 
variations  in  qSV-wave  amplitudes  (Figure  12).  In  the  anisotropic  medium,  there  is  a  weaker 
dipole  acting  in  the  vertical  direction,  but  at  the  same  time,  the  material  is  “softer”  in 
this  direction  (compare  Cu  and  C33,  Table  1).  Hence,  it  yields  more  easily  to  the  applied 
source,  and  the  overall  response  of  the  medium  is  nearly  isotropic,  at  least  compared  to  the 
isotropic  hard  rock  medium,  which  displays  the  effects  of  the  unequal  dipoles  in  the  source 
representation. 
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Transversely  Isotropic  Shale  Formation 

Another  type  of  anisotropic  medium  that  might  occur  in  a  nuclear  test  site  is  the  transverse 
isotropy  frequently  observed  in  shale  formations.  Since  the  axis  of  symmetry  in  these  forma¬ 
tions  is  generally  vertical,  our  simulation  of  the  radiation  from  a  finite  length  source  cavity 
would  correspond  to  a  model  of  a  vertical  shaft  in  the  shale  formation.  An  example  of  such 
a  formation  is  the  Pierre  Shale,  in  which  White  et  al.  (1983)  performed  a  VSP  experiment  to 
estimate  the  five  elastic  constants  of  the  medium.  One  of  their  estimates  is  shown  in  Table  2, 
and  the  phase  velocities  corresponding  to  these  elastic  constants  are  shown  in  Figure  14.  Un¬ 
like  the  hypothetical  fractured  medium  considered  above,  the  quasi-shear  wave  polarized  in 
the  plane  containing  the  cavity  axis  is  faster  than  the  horizontally  polarized  SH  mode  for 
angles  less  than  approximately  30  degrees.  In  other  words,  for  slowness  vector  angles  within 
30  degrees  of  the  cylinder  axis,  the  observed  qSV-wave  will  correspond  to  the  faster  of  the 
two  shear  wave  velocities  (qS2).  For  larger  angles,  the  qSV-wave  will  be  associated  with  the 
slower  phase  velocity  (qSi). 

Radiation  patterns  for  the  cavity  source  in  the  Pierre  Shale  display  some  of  the  same 
trends  as  did  the  source  in  the  fractured  medium  (Figure  15).  In  particular,  the  S-wave 
amplitude  increases  with  decreasing  cavity  aspect  ratio.  Like  the  cavities  in  the  fractured 
medium,  the  moment  tensor  representations  of  the  sources  (equation  (7))  are  almost  iden¬ 
tical  in  the  isotropic  low  velocity  formation  and  the  anisotropic  shale.  A  comparison  of  the 
patterns  for  the  isotropic  and  anisotropic  shales  also  shows  that  the  P-wave  radiation  pat¬ 
tern  is  slightly  less  isotropic  for  the  anisotropic  shale,  with  a  somewhat  diamond-like  shape 
(Figure  16).  The  most  significant  difference  from  the  patterns  computed  for  the  material 
with  aligned  fractures  (Figure  12)  is  that  the  qSV-wave  radiation  pattern  is  larger  relative  to 
the  qP-wave  pattern.  Another,  more  subtle,  feature  of  the  radiation  patterns  for  the  Pierre 
Shale  is  that  the  qSV-wave  radiation  patterns  for  aspect  ratio  5  show  small  lobes  oriented 
in  the  near  vertical  direction  (Figure  15). 
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Discussion 


The  results  presented  above  for  the  isotropic  medium  are,  in  general,  similar  to  far-field 
displacement  fields  obtained  in  earlier  work.  Glenn  et  al.  (1985)  used  finite  element  methods 
to  solve  for  displacement  fields  on  an  elhpsoidal  cavity  surface  and  applied  the  representation 
theorem  to  determine  the  far-field  displacements.  Subsequent  work  continued  this  effort  with 
an  analytic  model,  finding  that  when  aspect  ratio  is  10,  the  cylindrical  and  ellipsoidal  cavities 
yield  essentially  the  same  radiation  patterns  (Glenn  et  al.,  1986).  Likewise,  Rial  and  Moran 
(1986)  developed  an  analytical  formulation  representing  the  source  as  a  set  of  dipoles  (i.e.,  a 
moment  tensor  of  the  form  shown  in  equation  (7)).  The  relative  amplitudes  were  described  as 
a  function  of  aspect  ratio,  so  that  as  the  aspect  ratio  approached  one,  all  dipoles  had  the  same 
magnitude  (D„  =  Dfi)  to  yield  the  classical  result  for  the  spherical  cavity.  However,  Glenn 
et  al.(1986)  and  Rial  and  Moran  (1986)  both  found  for  large  aspect  ratios  that  =  2/ZDh. 
As  discussed  above,  the  BEM  results  for  aspect  ratio  5  correspond  to  =  9.9Dh,  a  much 
weaker  difference.  The  reason  for  this  disagreement  is  not  clear.  It  should  also  be  noted  that 
like  Rial  and  Moran  (1986),  we  find  that  the  radiation  patterns  do  not  change  significantly 
for  aspect  ratios  over  5. 

A  result  more  similar  to  ours  was  obtained  by  Stevens  et  al.  (1991),  who  simulated 
wavefields  radiated  from  ellipsoidal  cavities  of  aspect  ratio  4  by  coupling  a  non-linear,  near¬ 
field  finite  difference  computation  with  the  representation  theorem.  Applying  a  pressure  step 
to  the  cavity  wall,  they  found  a  radiation  pattern  that  showed  strong  variations  in  P-wave 
strength  when  frequencies  up  to  100  Hz  were  considered.  However,  filtering  the  signals  with 
an  anelastic  attenuation  operator  to  replicate  results  more  relevant  to  regional  observations 
yielded  signals  with  almost  no  variation  in  amplitude  with  direction.  This  behavior  is  very 
similar  to  our  results  obtained  with  the  BEM  for  cylindrical  cavities. 

Non-linear  effects,  which  are  neglected  in  the  BEM  computation,  may  also  strongly  in¬ 
fluence  the  far-fleld  radiation  pattern  of  an  explosion  source.  Stevens  et  al.  (1991)  showed 
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that  the  non-linear  interactions  of  an  air  shock  with  the  medium  surrounding  the  cavity  gen¬ 
erated  a  P-wave  radiation  pattern  that  was  even  more  isotropic  than  the  pattern  generated 
by  the  pressure  step  simulation.  This  suggests  that  the  BEM  results  for  the  instantaneous 
pressurization  model  of  the  explosion  should  be  considered  an  upper  bound  for  the  deviation 
of  the  radiation  pattern  from  that  of  the  isotropic,  spherical  source. 

The  neglect  of  non-linear  effects  in  the  BEM  simulation  could  be  worse  when  we  consider 
the  point  source  explosion  model,  especially  when  the  source  is  near  the  end  of  the  cavity 
(Figure  5).  It  is  difficult  to  predict  what  the  effect  of  such  non-linearities  on  the  radiation 
patterns  might  be,  but  the  results  discussed  above  for  the  sources  centered  in  the  cavity 
suggest  that  effect  might  be  to  reduce  signal  strength  variations.  Some  data  is  available 
which  helps  to  support  this  conclusion.  Although  field  data  where  source  cavity  shape  (for 
non-spherical  cavities),  source  location  and  waveforms  are  all  available  are  rather  scarce,  data 
from  Russian  experiments  in  the  1960’s  has  recently  become  available  (Murphy  et  al,  1995). 
Several  decoupling  tests  were  conducted  in  a  mine  in  the  Tywya  Mountains  of  Kirghizia, 
and  the  sources  were  located  in  both  spherical  and  elongated  cavities.  The  charges,  chemical 
explosives,  were  also  located  at  various  points  within  the  source  cavities.  Yields  of  the 
explosions  were  0.1,  1  and  6  tons.  Analysis  of  the  data  shows  that  there  is  little  variation  in 
peak  displacement  fields  with  changing  cavity  shape  or  source  location  (Murphy  et  al.,  1995). 
This  suggests  that  numerical  results  predicting  a  relatively  small  dependence  on  cylindrical 
cavity  aspect  ratio  are  good. 


Conclusions 

The  BEM  is  a  useful  approach  for  simulating  the  far-field  radiation  of  explosion  sources 
located  in  cylindrical  cavities.  We  have  applied  it  to  the  computation  of  radiation  patterns 
for  two  different  explosion  models,  a  volume  injection  point  source  and  an  instantaneous, 
uniform  pressurization  of  the  cavity.  While  the  wavefield  generated  by  the  point  source 
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strikes  the  cavity  wall  at  almost  constant  time,  the  amplitude  of  the  incident  field  is  inversely 
proportional  to  the  distance  from  the  source.  Our  computations  show  that  the  resulting 
radiation  patterns  are  sensitive  to  both  the  cavity  aspect  ratio  and  to  the  position  of  the 
source  along  the  cavity  axis.  When  the  source  is  near  the  end  of  the  cavity,  that  reradiation 
of  energy  from  this  near  end  dominates  the  far-field  signature  and  the  resulting  radiation 
pattern  resembles  that  of  a  point  force. 

On  the  other  hand,  the  instantaneous,  uniform  pressurization  model  of  the  source  gen¬ 
erates  radiation  patterns  that  are  similar  to  the  calculation  results  of  Stevens  et  al.  (1991) 
in  that  high  frequency  simulations  yield  P-wave  radiation  patterns  with  strong  variations  in 
amplitude.  The  corner  frequency  of  the  spectra  changes  with  direction  of  observation.  How¬ 
ever,  results  filtered  to  frequencies  more  typical  of  regional  wave  propagation  (1  Hz)  have 
lower  S-wave  amplitudes  and  P-wave  patterns  that  are  nearly  isotropic.  These  simulations 
suggest  that  the  location  of  a  source  in  a  tunnel-like  cavity  would  not  have  an  extremely 
strong  effect  on  regional  seismograms.  Application  of  this  model  to  two  transversely  isotropic 
media  shows  that  the  predicted  radiation  patterns  are  fairly  similar  to  the  isotropic  case. 
The  principal  difference  is  that  the  qSV-wave  amplitude  is  smallest  for  large  aspect  ratios. 
This  result  shows  that  anisotropy  does  not  necessarily  increase  the  strength  of  shear  wave 
energy  created  by  an  explosion  source. 
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Parameter 

Background  Value  (GPa) 

Perturbed  Value  (GPa) 

C'li 

35.2  GPa 

33.8  GPa 

C-ss 

35.2 

26.0 

C-ia 

13.9 

10.3 

<^44 

10.6 

9.53 

C-ee 

10.6 

10.6 

P 

2.2  g/cm® 

2.2 

Table  1:  Background  elastic  constants  and  perturbed  values  computed  using  the  theory  of 
Hudson  (1980,  1981).  Density  is  also  shown.  The  background  parameters  correspond  to  P 
and  S-wave  velocities  of  4000  m/s  and  2200  m/s,  respectively,  and  the  perturbed  values  were 
obtained  using  a  crack  density  of  0.05. 


Parameter 

Value 

Gii 

11.8  GPa 

Czz 

9.68 

Ciz 

7.09 

Cu 

1.62 

Gee 

2.16 

P 

2.25  g/cm^ 

Table  2:  Elastic  constants  and  density  of  Pierre  Shale  (White  et  al,  1983). 
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Figure  1:  Schematic  illustration  of  model  discretization  for  the  BEM.  In  general,  both  the 
ends  and  the  cylindrical  portion  of  the  cavity  are  discretized  as  shown  here.  An  example 
of  an  annulus  element  on  the  cavity  end  and  a  ring  element  on  the  cylinder  are  shaded 
in  the  figure. 
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Figure  2:  Synthetic  seismograms  computed  for  a  cylindrical  cavity  in  a  hard  rock  isotropic 
formation  with  P-wave  velocity  4000  m/s  and  S-wave  velocity  2200  m/s.  Other  param¬ 
eters  are  described  in  the  text.  Results  obtained  with  the  finite  difference  algorithm  are 
shown  in  the  left  hand  column,  while  BEM  results  are  on  the  right.  The  observation 
angle  for  these  seismograms  ranges  from  90  degrees  (perpendicular  to  the  cavity)  to  180 
degrees  (along  the  cavity  axis).  P-wave  signals  are  seen  on  the  radial  component,  while 
the  S-wave  appears  on  the  tangential  component. 


Figure  3:  Comparison  of  radiation  patterns  computed  using  finite  differences  and  the  BEM 
(Figure  2).  Cavity  orientation  in  this  figure  is  vertical. 
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Figure  4:  Radiation  patterns  created  by  a  volume  injection  point  source  at  the  center  of 
cavities  of  varying  aspect  ratio  located  in  the  hard  rock  isotropic  formation  (velocities  in 
Figure  2).  Since  the  dominant  source  frequency  was  specified  to  be  1  Hz,  these  patterns 
are  representations  of  energy  variations  that  might  be  observed  at  regional  distances. 


62 


2.0e-14 


40  M 


20M 


O.Oe-fOO 


-2.0e-14 


•2.0e-14 


O.Oe+00 


2.0e-14 


10  M 


1M 


Figure  5;  Radiation  patterns  created  by  a  volume  injection  point  source  at  varying  positions 
along  the  axis  of  an  80  m  long  cavity  with  aspect  ratio  5.  The  medium  is  the  hard  rock 
formation  (medium  properties  are  specified  in  Figure  2).  The  numbers  at  the  top  of  each 
plot  show  the  distance  of  the  source  from  the  end  of  the  cavity.  Therefore,  the  pattern 
at  the  upper  left  corresponds  to  the  source  at  the  center  of  the  cavity  and  is  identical  to 
the  plot  in  the  upper  left  corner  of  Figure  4. 
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Figure  6:  Log-log  plot  of  the  spectra  of  the  wavefields  generated  by  the  (A)  mstantaneous 
pressurization  and  (B)  volume  injection  point  source  models  in  the  cavity  of  aspect  ratio 
5  (isotropic  formation).  The  volume  injection  source  was  centered  in  the  cavity.  The  hard 
rock  formation  velocities  are  given  in  Figure  2.  Each  plot  shows  the  spectrum  observed 
at  3,  45  and  87  degrees  from  the  cavity  axis. 
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Figure  7:  High  frequency  radiation  patterns  of  the  (A)  instantaneous  pressurization  and  (B) 
point  source  sources  in  a  cavity  of  aspect  ratio  5  located  inside  the  hard  rock  isotropic 
formation  (length  40  m,  radius  4  m).  The  cavity  is  oriented  vertically  in  the  figure.  This 
figure  shows  the  amplitude  variation  in  synthetic  seismograms  computed  using  a  source 
having  a  step  function  time  dependence,  including  frequencies  up  to  100  Hz. 
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Figure  8:  Low  frequency  radiation  patterns  of  the  (A)  instantaneous  pressurization  and  (B) 
point  source  sources  in  a  cavity  of  aspect  ratio  5  located  inside  the  hard  rock  isotropic 
formation  (length  40  m,  radius  4  m).  The  cavity  is  oriented  vertically  in  the  figure.  This 
figure  shows  the  amplitude  variation  in  synthetic  seismograms  computed  using  a  Ricker 
wavelet  with  a  dominant  frequency  of  1  Hz. 
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Figure  9:  Radiation  patterns  for  the  instantaneous  pressurization  model  in  cavities  located 
in  the  hard  rock  isotropic  formation  (velocities  in  Figure  2).  In  all  cases,  the  cavity  length 
was  40  m,  and  it  is  oriented  vertically  in  the  figures.  The  frequency  domain  results  were 
convolved  with  a  Ricker  wavelet  of  1  Hz  to  simulate  the  frequency  ranges  seen  at  regional 
distances. 
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Figure  10:  Radiation  patterns  for  the  low  velocity  isotropic  shale.  The  figure  format  and 
cavity  dimensions  are  as  in  Figure  9.  Note  that  the  range  of  values  on  the  axes  has  been 
chosen  to  facifitate  direction  comparison  with  results  for  the  analogous  anisotropic  shale 
(Figure  15).  P-wave  velocity  of  this  soft  rock  formation  is  2290  m/s,  S-wave  velocity  is 
980  m/s  and  density  is  2250  kg/m^. 


68 


30 


45 


60 


75 


90 


0  15 

SLOWNESS  VECTOR  DIRECTION  (DEGREES) 


Figure  11:  Phase  velocities  of  the  fractured,  transversely  isotropic  medium  with  elastic 
constants  in  Table  1.  The  angle  of  the  slowness  vector  is  measured  from  the  2:-axis,  the 
axis  of  symmetry  of  the  medium. 
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Figure  12:  Radiation  patterns  for  the  instantaneous  pressurization  source  model  located  in 
the  fractured,  transversely  isotropic  formation.  Aspect  ratio  of  the  cavity  was  varied  as 
indicated  in  the  figure.  Medium  parameters  for  the  material  surrounding  the  cavity  are 
given  in  Table  1.  Note  that  the  fractures  are  aligned  perpendicular  to  the  cavity.  Because 
the  background  elastic  constants  used  to  compute  the  elastic  constants  for  the  medium 
were  those  of  the  isotropic  formation  considered  above,  a  comparison  of  these  radiation 
patterns  with  those  in  Figure  9  shows  the  influence  of  aligned  fractures  on  the  generation 
of  quasi-compressional  and  quasi-shear  wave  energy  by  the  explosion  source. 
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Figure  13:  Comparison  of  radiation  patterns  in  isotropic  hard  rock  formation  (left)  and 
anisotropic,  fractured  equivalent  (right). 
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Figure  14:  Phase  velocity  variation  in  the  transversely  isotropic  Pierre  shale.  Elastic  con¬ 
stants  were  estimated  by  White  et  al.  (1983)  using  VSP  data  (see  Table  2).  The  angle 
of  the  slowness  vector  is  measured  from  the  z— axis,  the  axis  of  symmetry. 
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Figure  15:  Radiation  patterns  for  the  instantaneous  pressurization  source  model  located  in 
cavities  insie  the  Pierre  Shale.  Aspect  ratio  of  the  40  m  long  cavity  varied  as  shown  in 
the  figure.  Medium  properties  are  specified  in  Table  2. 
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ASPECT  RATIO  1.25 


Figure  16:  Comparison  of  radiation  patterns  in  isotropic  low  velocity  formation  formation 
(left)  and  anisotropic  shale  (right).  The  velocities  in  the  isotropic  formation  (P-wave 
velocity  2290  m/s,  S-wave  velocity  980  m/s)  are  of  the  same  general  magnitude  as  those 
in  the  anisotropic  shale  (see  Figure  14). 
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